library(knitr)
library(rmdformats)
## Global options
options(max.print="75")
opts_chunk$set(#echo=FALSE,
cache=TRUE,
prompt=FALSE,
tidy=TRUE,
comment=NA,
message=FALSE,
warning=FALSE,
results='asis'
)
opts_knit$set(width=75)
print("notebook options imported successfully!")## [1] "notebook options imported successfully!"
## Load Required Libraries & Set Paths
library(tidyverse)
library(tidyquant)
library(lubridate)
library(timetk)
library(DT)
library(PerformanceAnalytics)
library(ROI.plugin.quadprog)
library(ROI.plugin.glpk)
# library(ROI.plugin.symphony)
setwd("~/Documents/GitHub/rob_gordon_2018")stocks <- c("SPY", "IYR", "QQQ", "TLT", "GLD", "IWM", "EEM", "DBC", "EFA")
stock_data <- tq_get(stocks
, get = "stock.prices"
, from = Sys.Date() - months(12)
, to = Sys.Date())
stock_data %>% head(10)| symbol | date | open | high | low | close | volume | adjusted |
|---|---|---|---|---|---|---|---|
| SPY | 2017-11-15 | 256.62 | 257.22 | 255.63 | 256.44 | 80811500 | 251.8421 |
| SPY | 2017-11-16 | 257.52 | 259.04 | 257.47 | 258.62 | 67777000 | 253.9830 |
| SPY | 2017-11-17 | 258.22 | 258.59 | 257.77 | 257.86 | 75756800 | 253.2366 |
| SPY | 2017-11-20 | 258.14 | 258.52 | 257.86 | 258.30 | 48075500 | 253.6688 |
| SPY | 2017-11-21 | 259.18 | 260.20 | 258.26 | 259.99 | 69176800 | 255.3284 |
| SPY | 2017-11-22 | 260.00 | 260.15 | 259.57 | 259.76 | 45033400 | 255.1026 |
| SPY | 2017-11-24 | 260.32 | 260.48 | 260.16 | 260.36 | 27856500 | 255.6918 |
| SPY | 2017-11-27 | 260.41 | 260.75 | 260.00 | 260.23 | 52274900 | 255.5642 |
| SPY | 2017-11-28 | 260.76 | 262.90 | 260.65 | 262.87 | 98971700 | 258.1568 |
| SPY | 2017-11-29 | 263.02 | 263.63 | 262.20 | 262.71 | 77512100 | 257.9997 |
# using ggplot2
stock_data %>%
filter(symbol %in% c("SPY", "GLD", "QQQ", "DBC")) %>%
ggplot(aes(x = date, y = adjusted, color = symbol)) +
geom_line(size = 1) +
geom_bbands(aes(high = high, low = low, close = close), ma_fun = SMA, n = 5,show.legend=TRUE) +
labs(title = "Daily Stock Prices",
x = "", y = "Adjusted Prices", color = "") +
facet_wrap(~ symbol, ncol = 2, scales = "free") +
scale_y_continuous(labels = scales::dollar) +
theme_tq() +
scale_color_tq()mo_returns <- stock_data %>%
group_by(symbol) %>%
tq_transmute(select = adjusted,
mutate_fun = periodReturn,
period = "monthly",
col_rename = "returns")
mo_returns %>%
mutate(returns = round(returns, 6)) %>%
arrange(desc(symbol), desc(date)) %>% datatable(class="compact", options = list(scrollX = TRUE, dom = "tp"))# simple returns = Rt = [ Pt - Pt-1 ]/ Pt-1
# simple_returns <-
# stock_data %>%
# group_by(symbol, year = year(date)) %>%
# summarise(adjusted = last(adjusted, order_by = year)) %>%
# mutate(return = (adjusted - lag(adjusted, order_by = year)) / lag(adjusted, order_by = year)
# , return = round(return, 6)) %>%
# arrange(desc(symbol), desc(year))
# # %>% filter(symbol == "SPY")
# datatable(simple_returns, class="compact", options = list(scrollX = TRUE, dom = "tp"))mo_returns %>%
filter(symbol %in% c("SPY", "GLD", "QQQ", "DBC")) %>%
ggplot(aes(x = date, y = returns, fill = symbol)) +
geom_line(aes(color=symbol)) +
geom_hline(yintercept = 0, color = palette_light()[[1]]) +
scale_y_continuous(labels = scales::percent) +
labs(title = " Monthly Returns",
subtitle = "stocks limited to fit",
y = "Monthly Returns", x = "") +
facet_wrap(~ symbol, ncol = 2,scales = "free") +
theme_tq() +
scale_fill_tq()mo_returns %>%
filter(symbol %in% c("SPY", "GLD", "QQQ", "DBC")) %>%
ggplot(aes(x = date, y = returns, fill = symbol)) +
geom_smooth(aes(color=symbol),method = 'loess' , formula = y ~ x, se=FALSE) +
geom_hline(yintercept = 0, color = palette_light()[[1]]) +
scale_y_continuous(labels = scales::percent) +
labs(title = " Monthly Returns",
subtitle = "stocks limited to fit",
y = "Monthly Returns", x = "") +
facet_wrap(~ symbol, ncol = 2,scales = "free") +
theme_tq() +
scale_fill_tq()init.investment <- 1000
growth <- mo_returns %>% arrange(date) %>%
mutate(final_value = init.investment * cumprod(1 + returns)) %>%
arrange(desc(final_value))
growth %>% filter(date == max(date)) %>% select(-date)| symbol | returns | final_value |
|---|---|---|
| QQQ | -0.0272053 | 1090.9494 |
| SPY | -0.0015889 | 1072.8945 |
| IWM | -0.0038659 | 1038.9866 |
| IYR | 0.0266291 | 1016.4514 |
| DBC | -0.0536557 | 996.2755 |
| GLD | -0.0044290 | 944.2385 |
| EFA | 0.0036824 | 941.7529 |
| TLT | 0.0069282 | 925.4958 |
| EEM | 0.0196629 | 893.9512 |
growth %>% ggplot(aes(x = date, y = final_value, color = symbol)) +
geom_line() +
# geom_smooth(method = "loess") +
labs(title = "Individual Portfolio: Comparing the Growth of $1K",
subtitle = "Quickly visualize performance",
x = "", y = "Investment Value") +
theme_tq() + theme(legend.position = "right") +
scale_y_continuous(labels = scales::dollar) ## winners (top 5)
winners <- growth %>% ungroup() %>% filter(date == max(date)) %>% mutate(rank = row_number()) %>%
top_n(5, final_value) %>% select(rank, symbol, final_value)
winners %>% arrange(desc(final_value)) %>% ggplot(aes(x = symbol, y = final_value,
fill = symbol)) + geom_bar(stat = "identity") + # geom_smooth(method = 'loess') +
labs(title = "Individual Portfolio: Comparing the Growth of $1K", subtitle = "Quickly visualize performance",
x = "", y = "Investment Value") + theme_tq() + theme(legend.position = "right") +
scale_y_continuous(labels = scales::dollar) + geom_text(aes(label = round(final_value,
2), vjust = 3, size = 10))winners$symbol[1] “QQQ” “SPY” “IWM” “IYR” “DBC”
# convert to ts object
stock_returns <- mo_returns %>% filter(symbol %in% winners$symbol) %>% spread(key = symbol,
value = returns) %>% tk_xts(date_var = date) %>% na.omit()
# stock_returns %>% class()hist_returns <- tq_get(winners$symbol, get = "stock.prices", from = Sys.Date() -
months(36), to = Sys.Date()) %>% group_by(symbol) %>% tq_transmute(select = adjusted,
mutate_fun = periodReturn, period = "daily", type = "arithmetic", col_rename = "returns")
hist_returns_ts <- hist_returns %>% spread(key = symbol, value = returns) %>%
tk_xts(date_var = date) %>% na.omit()# Markowitz Optimization involves minimizing the risk(variance) of returns and maximizing returns.
# https://rdrr.io/cran/PortfolioAnalytics/f/inst/doc/ROI_vignette.pdf
library(PortfolioAnalytics)
# Loading all suggested packages to fix ROI issue
# list.of.packages <- c("foreach", "DEoptim", "iterators", "fGarch", "Rglpk", "quadprog", "ROI", "ROI.plugin.glpk", "ROI.plugin.quadprog", "ROI.plugin.symphony", "pso", "GenSA", "corpcor", "testthat", "nloptr", "MASS", "robustbase")
# new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
# if(length(new.packages)) install.packages(new.packages)
#' Set up initial portfolio with basic constraints.
init.portf <- portfolio.spec(assets = colnames(hist_returns_ts))
# Add full investment constraint to the portfolio object
init.portf <- add.constraint(portfolio=init.portf, type="full_investment")
# Add long only constraints
init.portf <- add.constraint(portfolio=init.portf, type="box", min_sum=0.99, max_sum=1.01)
# Add objective to minimize variance
portf_minvar <- add.objective(portfolio=init.portf, type="risk", name="var")
# Dummy objectives for plotting and/or further analysis
portf_minvar <- add.objective(portf_minvar, type="return", name="mean", multiplier=0)
portf_minvar <- add.objective(portf_minvar, type="risk", name="ES", multiplier=0)
portf_minvar <- add.objective(portf_minvar, type="risk", name="StdDev", multiplier=0)
opt_gmv <- optimize.portfolio(R=hist_returns_ts,
portfolio=portf_minvar,
optimize_method="ROI",
trace=TRUE,
search_size=5000)
### Print Min Variance Achieved via Optimization ###
extractStats(opt_gmv)[1] ES
0.01513927
### Print Weights ###
weights <- extractWeights(opt_gmv)
weights_tbl <- as.tibble(as.list(weights)) %>% gather("symbol", "weight", 1:length(weights)) %>% arrange(desc(weight))
weights_tbl| symbol | weight |
|---|---|
| IYR | 0.4202591 |
| DBC | 0.3657620 |
| SPY | 0.2139789 |
| IWM | 0.0000000 |
| QQQ | 0.0000000 |
### Bar Chart Weights ###
opt_gmv %>% chart.Weights(plot.type="barplot", ylim=c(0,1), las=1, colorset = rainbow12equal, main="Optimal Weights", cex.axis = .8 )### Pie Chart Weights ###
library(scales)
pie <- ggplot(weights_tbl, aes(x="", y=weight, fill=symbol)) +
geom_bar(width = 1, stat = "identity") +
coord_polar("y", start=0) +
scale_fill_brewer("Blues") + #blank_theme() +
theme(axis.text.x=element_blank())+
geom_text(aes(y = weight/3 + c(0, cumsum(weight)[-length(weight)]),
label = percent(weight)), size=5)
piepaste("Min Variance Achieved: ", round(opt_gmv$objective_measures$ES,5))[1] “Min Variance Achieved: 0.01514”
### Chart Risk Reward ###
chart.RiskReward(opt_gmv, rp=TRUE, risk.col = "ES", chart.assets=TRUE)annualized.moments <- function(R, scale=12, portfolio=NULL){
out <- list()
out$mu <- matrix(colMeans(R), ncol=1)
out$sigma <- cov(R)/scale
return(out)
}
### CALCULATE & PLOT EFFICIENT FRONTIER ###
prt_ef <- create.EfficientFrontier(R=12*hist_returns_ts, portfolio=portf_minvar, type="mean-StdDev",
match.col = "StdDev", momentFUN="annualized.moments", scale=12)
### Chart Efficient Frontier ###
xlim <- range(prt_ef$frontier[,2])*c(1, 1.5)
ylim <- range(prt_ef$frontier[,1])*c(.80, 1.05)
chart.EfficientFrontier(prt_ef, match.col="StdDev", chart.assets = TRUE, n.portfolios = 1000
, labels.assets = TRUE, xlim=xlim, ylim=ylim, main = "Efficient Frontier")
points(with(annualized.moments(12*stock_returns, scale=12), cbind(sqrt(diag(sigma)), mu)), pch=19 )
text(with(annualized.moments(12*stock_returns, scale=12), cbind(sqrt(diag(sigma)), mu)),
labels=colnames(stock_returns), cex=.8, pos=4) ### Chart EF Weights ###
chart.EF.Weights(prt_ef, match.col="StdDev", colorset = rainbow12equal, main = "Efficient Frontier")# https://rdrr.io/cran/PortfolioAnalytics/f/inst/doc/ROI_vignette.pdf
bt_gmv <- optimize.portfolio.rebalancing(R=hist_returns_ts, portfolio=portf_minvar,
optimize_method="ROI",
trace = TRUE,
rebalance_on="months",
training_period=36,
# trailing_period=trailing,
rolling_window = 12
)
summary(bt_gmv)$annualized_returns[[1]][1] 0.1133622
summary(bt_gmv)$annualized_StdDev[[1]][1] 0.1018277
bt_weights <- summary(bt_gmv)$weights
bt_ret <- summary(bt_gmv)$portfolio_returns
charts.PerformanceSummary(cbind(bt_ret, opt_gmv$R), main="Optimization Performance", method = "StdDev", wealth.index = TRUE, begin = "axis")x <- Return.portfolio(hist_returns_ts, bt_weights, wealth.index = TRUE, verbose = TRUE)
chart.CumReturns(x$returns, col=bluemono, begin="axis")chart.StackedBar(x$BOP.Weight)chart.StackedBar(x$BOP.Value)# CHART BACKTEST WEIGHTS
bt_gmv %>% chart.Weights(ylim=c(0,1), las=1, main="Optimal Weights", cex.axis = 3 ,cex.legend = 3, col=rainbow12equal) # , legend.loc = "topright"# ret.bt.gmv <- do.call(cbind, lapply(bt_gmv, function(x) summary(bt_gmv)$portfolio_returns))
# colnames(ret.bt.gmv) <- c("min ES", "min ES RB", "min ES Eq RB", "stuff", "other stuff")
Rb <- Return.portfolio(hist_returns_ts, weights = c(.2,.2,.2,.2,.2), wealth.index = TRUE, verbose = TRUE)
chart.RelativePerformance(as.xts(x$returns), as.xts(Rb$returns), colorset = rainbow12equal, main = "Performance Relative to Even Dist")| symbol | ArithmeticMean | GeometricMean | Kurtosis | LCLMean(0.95) | Maximum | Median | Minimum | NAs | Observations | Quartile1 | Quartile3 | SEMean | Skewness | Stdev | UCLMean(0.95) | Variance |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| SPY | 0.0060 | 0.0054 | 0.0073 | -0.0146 | 0.0564 | 0.0059 | -0.0691 | 0 | 13 | -0.0016 | 0.0319 | 0.0094 | -0.7344 | 0.0340 | 0.0265 | 0.0012 |
| IYR | 0.0017 | 0.0013 | -0.4254 | -0.0176 | 0.0406 | 0.0023 | -0.0666 | 0 | 13 | -0.0239 | 0.0266 | 0.0089 | -0.6159 | 0.0319 | 0.0210 | 0.0010 |
| QQQ | 0.0077 | 0.0067 | -0.0764 | -0.0197 | 0.0876 | 0.0060 | -0.0860 | 0 | 13 | -0.0129 | 0.0280 | 0.0126 | -0.2076 | 0.0454 | 0.0351 | 0.0021 |
| TLT | -0.0057 | -0.0059 | -1.5196 | -0.0190 | 0.0286 | -0.0114 | -0.0326 | 0 | 13 | -0.0286 | 0.0131 | 0.0061 | 0.1430 | 0.0220 | 0.0076 | 0.0005 |
| GLD | -0.0042 | -0.0044 | -0.7278 | -0.0162 | 0.0323 | -0.0066 | -0.0361 | 0 | 13 | -0.0208 | 0.0063 | 0.0055 | 0.3733 | 0.0199 | 0.0078 | 0.0004 |
| IWM | 0.0039 | 0.0029 | 1.3954 | -0.0230 | 0.0616 | 0.0098 | -0.1099 | 0 | 13 | -0.0040 | 0.0256 | 0.0124 | -1.1468 | 0.0446 | 0.0308 | 0.0020 |
| EEM | -0.0076 | -0.0086 | -0.4608 | -0.0354 | 0.0830 | -0.0058 | -0.0876 | 0 | 13 | -0.0377 | 0.0197 | 0.0128 | 0.1761 | 0.0461 | 0.0202 | 0.0021 |
| DBC | 0.0002 | -0.0003 | -1.1508 | -0.0198 | 0.0342 | 0.0075 | -0.0562 | 0 | 13 | -0.0243 | 0.0278 | 0.0092 | -0.5503 | 0.0331 | 0.0202 | 0.0011 |
| EFA | -0.0040 | -0.0046 | 0.2145 | -0.0249 | 0.0502 | 0.0037 | -0.0813 | 0 | 13 | -0.0189 | 0.0152 | 0.0096 | -0.7049 | 0.0345 | 0.0168 | 0.0012 |
| symbol | AnnualizedReturn | AnnualizedSharpe(Rf=0%) | AnnualizedStdDev |
|---|---|---|---|
| SPY | 0.0671 | 0.5690 | 0.1179 |
| IYR | 0.0152 | 0.1371 | 0.1107 |
| QQQ | 0.0837 | 0.5325 | 0.1571 |
| TLT | -0.0690 | -0.9040 | 0.0763 |
| GLD | -0.0516 | -0.7486 | 0.0689 |
| IWM | 0.0359 | 0.2328 | 0.1543 |
| EEM | -0.0983 | -0.6160 | 0.1596 |
| DBC | -0.0034 | -0.0300 | 0.1146 |
| EFA | -0.0539 | -0.4514 | 0.1194 |
| symbol | DownsideDeviation(0%) | DownsideDeviation(MAR=10%) | DownsideDeviation(Rf=0%) | GainDeviation | HistoricalES(95%) | HistoricalVaR(95%) | LossDeviation | MaximumDrawdown | ModifiedES(95%) | ModifiedVaR(95%) | SemiDeviation |
|---|---|---|---|---|---|---|---|---|---|---|---|
| SPY | 0.0230 | 0.0269 | 0.0230 | 0.0178 | -0.0691 | -0.0495 | 0.0279 | 0.0706 | -0.0673 | -0.0543 | 0.0257 |
| IYR | 0.0228 | 0.0273 | 0.0228 | 0.0146 | -0.0666 | -0.0448 | 0.0235 | 0.0959 | -0.0637 | -0.0542 | 0.0237 |
| QQQ | 0.0277 | 0.0318 | 0.0277 | 0.0303 | -0.0860 | -0.0589 | 0.0324 | 0.1133 | -0.0861 | -0.0666 | 0.0315 |
| TLT | 0.0185 | 0.0244 | 0.0185 | 0.0085 | -0.0326 | -0.0313 | 0.0084 | 0.0868 | -0.0430 | -0.0403 | 0.0145 |
| GLD | 0.0152 | 0.0213 | 0.0152 | 0.0107 | -0.0361 | -0.0279 | 0.0109 | 0.1166 | -0.0382 | -0.0339 | 0.0124 |
| IWM | 0.0330 | 0.0367 | 0.0330 | 0.0217 | -0.1099 | -0.0670 | 0.0438 | 0.1339 | -0.1043 | -0.0782 | 0.0346 |
| EEM | 0.0353 | 0.0406 | 0.0353 | 0.0282 | -0.0876 | -0.0704 | 0.0263 | 0.2276 | -0.0919 | -0.0786 | 0.0306 |
| DBC | 0.0246 | 0.0295 | 0.0246 | 0.0117 | -0.0562 | -0.0547 | 0.0171 | 0.1098 | -0.0613 | -0.0576 | 0.0247 |
| EFA | 0.0279 | 0.0326 | 0.0279 | 0.0154 | -0.0813 | -0.0615 | 0.0275 | 0.1373 | -0.0803 | -0.0647 | 0.0258 |
| symbol | Annualiseddownsiderisk | Downsidepotential | Monthlydownsiderisk | Omega | Omega-sharperatio | Sortinoratio | Upsidepotential | Upsidepotentialratio |
|---|---|---|---|---|---|---|---|---|
| SPY | 0.0795 | 0.0103 | 0.0230 | 1.5770 | 0.5770 | 0.2600 | 0.0163 | 0.5693 |
| IYR | 0.0788 | 0.0116 | 0.0228 | 1.1497 | 0.1497 | 0.0761 | 0.0133 | 0.6734 |
| QQQ | 0.0959 | 0.0131 | 0.0277 | 1.5875 | 0.5875 | 0.2769 | 0.0207 | 0.7542 |
| TLT | 0.0640 | 0.0129 | 0.0185 | 0.5567 | -0.4433 | -0.3093 | 0.0072 | 0.6175 |
| GLD | 0.0526 | 0.0104 | 0.0152 | 0.5959 | -0.4041 | -0.2778 | 0.0062 | 1.1078 |
| IWM | 0.1142 | 0.0138 | 0.0330 | 1.2820 | 0.2820 | 0.1181 | 0.0177 | 0.5410 |
| EEM | 0.1221 | 0.0223 | 0.0353 | 0.6591 | -0.3409 | -0.2157 | 0.0147 | 0.6629 |
| DBC | 0.0852 | 0.0141 | 0.0246 | 1.0159 | 0.0159 | 0.0091 | 0.0143 | 0.5856 |
| EFA | 0.0967 | 0.0150 | 0.0279 | 0.7306 | -0.2694 | -0.1449 | 0.0110 | 0.4958 |
| symbol | VaR |
|---|---|
| SPY | -0.0543204 |
| IYR | -0.0541734 |
| QQQ | -0.0666155 |
| TLT | -0.0403022 |
| GLD | -0.0338613 |
| IWM | -0.0782136 |
| EEM | -0.0785758 |
| DBC | -0.0575836 |
| EFA | -0.0646908 |
require(quadprog)
# min_x(-d^T x + 1/2 b^T D x) r.t A.x>=b
MV_QP <- function(nx, tarRet, Sig = NULL, long_only = FALSE) {
if (is.null(Sig))
Sig = cov(nx)
dvec = rep(0, ncol(Sig))
meq = 2
Amat = rbind(rep(1, ncol(Sig)), apply(nx, 2, mean))
bvec = c(1, tarRet)
if (long_only) {
meq = 1
Amat = Amat[-1, ]
Amat = rbind(Amat, diag(1, ncol(Sig)), rep(1, ncol(Sig)), rep(-1, ncol(Sig)))
bvec = bvec[-1]
bvec = c(bvec, rep(0, ncol(Sig)), 0.98, -1.02)
}
sol <- solve.QP(Dmat = Sig, dvec, t(Amat), bvec, meq = meq)$solution
}
steps = 50
x = stock_returns
µ.b <- apply(X = x, 2, FUN = mean)
long_only = TRUE
range.bl <- seq(from = min(µ.b), to = max(µ.b) * ifelse(long_only, 1, 1.6),
length.out = steps)
risk.bl <- t(sapply(range.bl, function(targetReturn) {
w <- MV_QP(x, targetReturn, long_only = long_only)
c(sd(x %*% w), w)
}))
weigthsl = round(risk.bl[, -1], 4)
colnames(weigthsl) = colnames(x)
weigthsl %>% head() DBC IWM IYR QQQ SPY
[1,] 0.9770 0 0.0030 0 0.0000 [2,] 0.8763 0 0.1037 0 0.0000 [3,] 0.7756 0 0.2044 0 0.0000 [4,] 0.6749 0 0.3051 0 0.0000 [5,] 0.5741 0 0.4059 0 0.0000 [6,] 0.4908 0 0.4831 0 0.0062
risk.bl = risk.bl[, 1]
rets.bl = weigthsl %*% µ.b
fan = 12
plot(x = risk.bl * fan^0.5, y = rets.bl * fan, col = 2, pch = 21, xlab = "Annualized Risk ",
ylab = "Annualized Return", main = "long only EF with solve.QP")# demo(multi_layer_optimization)